library(qiime2R)
library(ggplot2)
library(vegan)
## Warning: package 'vegan' was built under R version 4.1.3
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.1.3
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.1.3
## This is vegan 2.6-4
library(plyr)
## Warning: package 'plyr' was built under R version 4.1.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(scales)
## Warning: package 'scales' was built under R version 4.1.3
library(grid)
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.1.3
library(phyloseq)
library(picante)
## Warning: package 'picante' was built under R version 4.1.3
## Loading required package: ape
## Warning: package 'ape' was built under R version 4.1.3
##
## Attaching package: 'ape'
## The following object is masked from 'package:dplyr':
##
## where
## Loading required package: nlme
## Warning: package 'nlme' was built under R version 4.1.3
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.3
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
##
## smiths
library(viridis)
## Loading required package: viridisLite
##
## Attaching package: 'viridis'
## The following object is masked from 'package:scales':
##
## viridis_pal
library(qiime2R)
library(DESeq2)
## Loading required package: S4Vectors
## Warning: package 'S4Vectors' was built under R version 4.1.3
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:tidyr':
##
## expand
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:plyr':
##
## rename
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:nlme':
##
## collapse
## The following object is masked from 'package:phyloseq':
##
## distance
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
## The following object is masked from 'package:plyr':
##
## desc
## The following object is masked from 'package:grDevices':
##
## windows
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## Warning: package 'matrixStats' was built under R version 4.1.3
##
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
##
## count
## The following object is masked from 'package:plyr':
##
## count
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
## The following object is masked from 'package:phyloseq':
##
## sampleNames
library(patchwork)
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 4.1.3
library(microViz)
## microViz version 0.10.10 - Copyright (C) 2023 David Barnett
## ! Website: https://david-barnett.github.io/microViz
## v Useful? For citation details, run: `citation("microViz")`
## x Silence? `suppressPackageStartupMessages(library(microViz))`
library(speedyseq)
##
## Attaching package: 'speedyseq'
## The following objects are masked from 'package:phyloseq':
##
## filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
## tip_glom, transform_sample_counts
library(ComplexHeatmap)
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(ggVennDiagram)
## Warning: package 'ggVennDiagram' was built under R version 4.1.3
library(SuperExactTest)
##
## Attaching package: 'SuperExactTest'
## The following objects are masked from 'package:GenomicRanges':
##
## intersect, union
## The following object is masked from 'package:GenomeInfoDb':
##
## intersect
## The following objects are masked from 'package:IRanges':
##
## intersect, union
## The following objects are masked from 'package:S4Vectors':
##
## intersect, union
## The following objects are masked from 'package:BiocGenerics':
##
## intersect, union
## The following objects are masked from 'package:dplyr':
##
## intersect, union
## The following objects are masked from 'package:base':
##
## intersect, union
library(nVennR)
library(vegan)
library(plyr)
library(dplyr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)
library(picante)
library(tidyr)
library(viridis)
library(qiime2R)
library(DESeq2)
library(patchwork)
library(RColorBrewer)
library(speedyseq)
library(RColorBrewer)
library(microViz)
library(ComplexHeatmap)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ComplexHeatmap':
##
## add_heatmap
## The following object is masked from 'package:microViz':
##
## add_paths
## The following object is masked from 'package:IRanges':
##
## slice
## The following object is masked from 'package:S4Vectors':
##
## rename
## The following objects are masked from 'package:plyr':
##
## arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(seecolor)
library(htmlwidgets)
## Warning: package 'htmlwidgets' was built under R version 4.1.3
library(htmltools)
## Warning: package 'htmltools' was built under R version 4.1.3
set.seed(17384)
PRphyseq <- qza_to_phyloseq(features="S003-S015-table.qza", tree="S003-S015-rooted-tree.qza",taxonomy="S003-S015-taxonomy.qza", metadata = "S003-S015-Field-Measurements-2.txt")
rank_names(PRphyseq)
## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"
sample_variables(PRphyseq)
## [1] "Project" "Sub_Location"
## [3] "Location" "Season"
## [5] "Month" "Year"
## [7] "Sample_Type" "Sample_ID"
## [9] "Day" "Station_ID"
## [11] "Station_Name" "Collected_by"
## [13] "Time" "Samples_Collected_by"
## [15] "Temperature_C" "mm_Hg"
## [17] "DO_mg_L" "DO_LDO"
## [19] "Specific_Conductance_uS_cm" "Salinity_ppt"
## [21] "pH" "pH_mv"
## [23] "IVCH_ug_L" "Turbidity_.Fluorometer"
## [25] "Comments"
## Removing chloroplast sequences and any contaminant sequences
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Eukaryota")
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Archaea")
PRphyseq <- subset_taxa(PRphyseq, Order != "o__Chloroplast")
PRphyseq <- subset_taxa(PRphyseq, Order != "Chloroplast")
PRphyseq <- subset_taxa(PRphyseq, Family != "f__Mitochondria")
PRphyseq <- subset_taxa(PRphyseq, Family != "Mitochondria")
PRphyseq <- subset_taxa(PRphyseq, Family != "Chloroplast")
##Check that contaminant sequences are removed (easiest to save as data frame and search) taxtabl <- as.data.frame(tax_table(S003physeq)) #At this point, blanks and positive controls should also be removed #S003physeq = subset_samples(S003physeq, Sample != "mock")
View((tax_table(PRphyseq)))
##This will allow you to open your data frame. You can search in the top right corner of the dataframe for the sequences you wanted to remove in above code to ensure they were removed.
PRphyseq <- PRphyseq %>% subset_samples(Project %in% c("S003","S015")) ##This will create a subset (based on project column) of the samples we would like to use in the analysis. If you would like to keep all samples in the metadata file, you can ignore this step.
PR_Normal <- PRphyseq %>% subset_samples(Sample_Type %in% c("N"))
readsumsdf = data.frame(nreads = sort(taxa_sums(PR_Normal), TRUE), sorted = 1:ntaxa(PR_Normal), type = "ASVs")
readsumsdf = rbind(readsumsdf, data.frame(nreads = sort(sample_sums(PR_Normal), TRUE), sorted = 1:nsamples(PR_Normal), type = "Samples"))
title = "Total number of reads"
p = ggplot(readsumsdf, aes(x = sorted, y = nreads)) + geom_bar(stat = "identity")
p + ggtitle(title) + scale_y_log10() + facet_wrap(~type, 1, scales = "free")
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 1908 rows containing missing values (`geom_bar()`).
Plot a rarefaction curve using vegan:
taxa_are_rows(PR_Normal)
## [1] TRUE
mat <- t(otu_table(PR_Normal))
class(mat) <- "matrix"
## Warning in class(mat) <- "matrix": Setting class(x) to "matrix" sets attribute
## to NULL; result will no longer be an S4 object
class(mat)
## [1] "matrix" "array"
mat <- as(t(otu_table(PR_Normal)), "matrix")
class(mat)
## [1] "matrix" "array"
raremax <- min(rowSums(mat))
system.time(rarecurve(mat, step = 100, sample = raremax, col = "blue", label = FALSE))
## user system elapsed
## 4.28 0.17 6.32
#Sum reads and sort by samples with least to greatest number of reads
samples_reads <- sort(sample_sums(PR_Normal))
#Create table
knitr::kable(samples_reads) %>%
kableExtra::kable_styling("striped", latex_options="scale_down") %>%
kableExtra::scroll_box(width = "100%")
| x | |
|---|---|
| PR109C | 14 |
| pr242a | 679 |
| PR325 | 699 |
| PR164C | 700 |
| PR151C | 855 |
| PR566A | 1040 |
| PR321 | 1088 |
| PR147C | 1128 |
| PR564A | 1219 |
| lsj3-1 | 1251 |
| PR323 | 1323 |
| PR559A | 1335 |
| PR145C | 1351 |
| PR422 | 1363 |
| PR425 | 1372 |
| PR162C | 1450 |
| PR141C | 1461 |
| PR421 | 1479 |
| PR317 | 1481 |
| PR399 | 1579 |
| PR561A | 1624 |
| PR311 | 1634 |
| PR420 | 1700 |
| PR139C | 1730 |
| PR333 | 1785 |
| PR558A | 1793 |
| PR554A | 1972 |
| PR419 | 2037 |
| PR417 | 2075 |
| PR160C | 2191 |
| PR107C | 2201 |
| PR556A | 2311 |
| PR334 | 2388 |
| PR557A | 2532 |
| PR392 | 2952 |
| PR395 | 3001 |
| PR149C | 3054 |
| PR394 | 3055 |
| PR332 | 3152 |
| PR393 | 3167 |
| PR555A | 3238 |
| PR416 | 3539 |
| PR562A | 3614 |
| PR415 | 3728 |
| PR328 | 3792 |
| PR400 | 3815 |
| PR316 | 4073 |
| PR565A | 4122 |
| PR331 | 4337 |
| PR402 | 4429 |
| PR560A | 4453 |
| PR401 | 4461 |
| PR563A | 4999 |
| PR315 | 5085 |
| 239 | 5356 |
| 18 | 5464 |
| pr248a | 5585 |
| PR318 | 5696 |
| PR102C | 6357 |
| PR322 | 6779 |
| pr250a | 7781 |
| pr262a | 8718 |
| pr271a | 9101 |
| 101 | 9859 |
| pr267a | 10277 |
| 220 | 10514 |
| pr247a | 10711 |
| 44 | 10755 |
| 104 | 10896 |
| cs1-1 | 11179 |
| cmpw-3 | 11765 |
| pr295a | 12285 |
| lsj3-2 | 12316 |
| cmpe-3 | 12325 |
| pr249a | 12448 |
| cmp2-2 | 12550 |
| 142 | 12638 |
| 198 | 12984 |
| lp1-2 | 13656 |
| 42 | 13896 |
| cmp2-1 | 13943 |
| cmpe-2 | 14167 |
| cmp2-3 | 14339 |
| lp1-3 | 14775 |
| 5 | 14776 |
| 218 | 15212 |
| 181 | 15472 |
| cmpe-1 | 15558 |
| cmpw-1 | 15630 |
| 183 | 16414 |
| 63 | 16771 |
| lsj3-3 | 16802 |
| cs1-2 | 19707 |
| 95 | 21854 |
| lp1-1 | 23460 |
| 106 | 25062 |
| 90 | 25753 |
| cs1-3 | 29344 |
| cmpw-2 | 51074 |
| pr293a | 73188 |
#Remove samples with less than 1,000 reads
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR105C")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR109C")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "pr242a")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR325")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR164C")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR151C")
#For each ASV, find the number of samples in which each ASV is 0, then divide by the total number of samples
test_filter_method <- as.data.frame(rowSums(otu_table(PR_Normal) == 0)/ncol(otu_table(PR_Normal)))
#Change the name of the column in the test_filter_method4 dataframe; this column contains the proportion of samples in which each ASV is 0
colnames(test_filter_method) <- c('samplew0')
#Create a bar plot
ggplot(data = test_filter_method) + geom_bar(mapping = aes(x= samplew0)) +labs(x = "Proportion of samples in which ASV = 0", y = "# of ASV's") +theme(text = element_text(size = 18), axis.title = element_text(size = 15),panel.spacing = unit(1, "lines"),panel.border = element_rect(colour = "black", fill = NA, size = 0.5), panel.background = element_blank())
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
##How many singletons and doubletons are present
#Create a dataframe that includes only the total number of times that each ASV occurs
readsumsdf_ASVs <- readsumsdf %>%
filter(type == "ASVs")
#How many singletons are present?
length(which(readsumsdf_ASVs$nreads <= 0))
## [1] 1908
length(which(readsumsdf_ASVs$nreads == 1))
## [1] 7
#How many doubletons are present?
length(which(readsumsdf_ASVs$nreads == 2))
## [1] 158
#Filter out low abundance ASVs
PR_Normal = prune_taxa(taxa_sums(PR_Normal) > 2, PR_Normal)
Examine the abundance of individual ASVs again:
#For each ASV, find the number of samples in which each ASV is 0, then divide by the total number of samples
test_filter_method_filtered <- as.data.frame(rowSums(otu_table(PR_Normal) == 0)/ncol(otu_table(PR_Normal)))
#Change the name of the column in the test_filter_method4 dataframe; this column contains the proportion of samples in which each ASV is 0
colnames(test_filter_method_filtered) <- c('samplew0')
#Create a bar plot
ggplot(data = test_filter_method_filtered) + geom_bar(mapping = aes(x= samplew0)) +
labs(x = "Proportion of samples in which ASV = 0", y = "# of ASV's") +
theme(text = element_text(size = 18), axis.title = element_text(size = 15),
panel.spacing = unit(1, "lines"),
panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_blank())
sample_data(PR_Normal)$Sub_Location <- factor(sample_data(PR_Normal)$Sub_Location, levels = c("BSJ1","BSJ1-ENR","BSJ2","BSJ2-ENR","BSJ3","BSJ3-ENR","BSJ4-ENR","BSJ5-ENR","BSJ6-ENR","BSJ7-ENR","CMP1","CMP1-ENR","CMP2","CMP2-ENR","CMPE","CMPW","CS1","CS1-ENR","CS2","CS2-ENR","CSA1-ENR","LLC1","LLC1-ENR","LP1","LP1-ENR","LP1-ENR","LP1L-ENR","LP2-ENR","LP3-ENR","LP4-ENR","LP5-ENR","LSJ1","LSJ1-ENR","LSJ2","LSJ2-ENR","LSJ3","LT1","LT1-ENR","LT2","LT2-ENR","LT3-ENR","QB1-ENR","QSA1-ENR","RPN1-ENR"),
labels = c("BSJ1","BSJ1-ENR","BSJ2","BSJ2-ENR","BSJ3","BSJ3-ENR","BSJ4-ENR","BSJ5-ENR","BSJ6-ENR","BSJ7-ENR","CMP1","CMP1-ENR","CMP2","CMP2-ENR","CMPE","CMPW","CS1","CS1-ENR","CS2","CS2-ENR","CSA1-ENR","LLC1","LLC1-ENR","LP1","LP1-ENR","LP1-ENR","LP1L-ENR","LP2-ENR","LP3-ENR","LP4-ENR","LP5-ENR","LSJ1","LSJ1-ENR","LSJ2","LSJ2-ENR","LSJ3","LT1","LT1-ENR","LT2","LT2-ENR","LT3-ENR","QB1-ENR","QSA1-ENR","RPN1-ENR"))
sample_data(PR_Normal)$Location <- factor(sample_data(PR_Normal)$Location, levels = c(
"BSJ","BSJ-ENR","CMP","CMP-ENR","CMPE","CMPW", "CS","CS-ENR", "CSA","CSA-ENR","LLC","LLC-ENR","LP","LP-ENR","LSJ","LSJ-ENR","LT","LT-ENR","QB","QB-ENR","QSA","QSA-ENR","RPN","RPN-ENR"),
labels = c("BSJ","BSJ-ENR","CMP","CMP-ENR","CMPE","CMPW", "CS","CS-ENR", "CSA","CSA-ENR","LLC","LLC-ENR","LP","LP-ENR","LSJ","LSJ-ENR","LT","LT-ENR","QB","QB-ENR","QSA","QSA-ENR","RPN","RPN-ENR"))
sample_data(PR_Normal)$Season <- factor(sample_data(PR_Normal)$Season, levels = c("W", "D"),
labels = c("W", "D"))
sample_data(PR_Normal)$Project <- factor(sample_data(PR_Normal)$Project,
levels = c("S003","S015"), labels = c("S003","S015"))
sample_data(PR_Normal)$Month <- factor(sample_data(PR_Normal)$Month,
levels = c("April","August","December","January","July","June","March","May","November","October","September"), labels = c("April","August","December","January","July","June","March","May","November","October","September"))
sample_data(PR_Normal)$Year <- factor(sample_data(PR_Normal)$Year,
levels = c("2021","2022","2023"), labels = c("2021","2022","2023"))
sample_data(PR_Normal)$Sample_Type <- factor(sample_data(PR_Normal)$Sample_Type,
levels = c("N","E"), labels = c("N","E"))
sample_data(PR_Normal)$Sample_ID <- factor(sample_data(PR_Normal)$Sample_ID,
levels = c("PR419","PR417","PR422","PR421","PR420","PR416","PR415","90","44","42","63","5","18","pr293a", "pr271a","pr262a","pr285a","pr295a","pr267a","pr290a","PR425","PR325","PR323","PR321","PR322","PR332","PR333","PR331","PR315","PR318","PR328","PR317","PR316","PR334", "PR311","PR563A","PR562A","PR560A","PR561A","PR565A","cmp2-3","cmpe-3","cmpw-3","cs1-3","PR556A","PR559A","PR564A","lp1-3","PR558A","PR557A","PR566A","lsj3-3","PR555A","PR554A","cmp2-2","cmpe-2","cmpw-2","cs1-2","lp1-2","lsj3-2","PR395","PR394","PR392", "PR393","PR400","PR401","PR399","PR402","PR464A","PR465A","PR466A","PR467A","PR468A","PR469A","PR470A","PR473A","cmp2-1","PR472A","cmpe-1","cmpw-1","cs1-1","PR476A","PR477A","PR474A","PR478A","lp1-1","PR488A","PR489A","PR490A","PR491A","PR492A","PR493A","PR479A","PR480A","lsj3-1","PR481A","PR482A","PR483A","PR484A","PR485A","PR471A","Mock","Undeter","B53","B84","B59","B43","B36","B44","B107","B90","B126","B141","B128","B152","B134","B132","B140","B89","B100","B102","B60","B76","B56","blank","blank2","Mock-com","B130","B136","B138","B33","B25","B41","B81","B68","B57","B87","B91","B117","B143","B139","B154","B69","B52","B75","B99","B118","B124","B26","B34","B42","pr249a","pr248a","pr250a","pr242a","pr247a","183","218","181","198","239","220","PR160C","PR162C","PR164C","PR102C","101","PR109C","108","PR107C","106","PR143C","142","PR145C","PR151C","PR96C","95","PR147C","PR149C","PR105C","104","PR139C","PR141C"
),
labels = c("PR419","PR417","PR422","PR421","PR420","PR416","PR415","90","44","42","63","5","18","pr293a", "pr271a","pr262a","pr285a","pr295a","pr267a","pr290a","PR425","PR325","PR323","PR321","PR322","PR332","PR333","PR331","PR315","PR318","PR328","PR317","PR316","PR334", "PR311","PR563A","PR562A","PR560A","PR561A","PR565A","cmp2-3","cmpe-3","cmpw-3","cs1-3","PR556A","PR559A","PR564A","lp1-3","PR558A","PR557A","PR566A","lsj3-3","PR555A","PR554A","cmp2-2","cmpe-2","cmpw-2","cs1-2","lp1-2","lsj3-2","PR395","PR394","PR392", "PR393","PR400","PR401","PR399","PR402","PR464A","PR465A","PR466A","PR467A","PR468A","PR469A","PR470A","PR473A","cmp2-1","PR472A","cmpe-1","cmpw-1","cs1-1","PR476A","PR477A","PR474A","PR478A","lp1-1","PR488A","PR489A","PR490A","PR491A","PR492A","PR493A","PR479A","PR480A","lsj3-1","PR481A","PR482A","PR483A","PR484A","PR485A","PR471A","Mock","Undeter","B53","B84","B59","B43","B36","B44","B107","B90","B126","B141","B128","B152","B134","B132","B140","B89","B100","B102","B60","B76","B56","blank","blank2","Mock-com","B130","B136","B138","B33","B25","B41","B81","B68","B57","B87","B91","B117","B143","B139","B154","B69","B52","B75","B99","B118","B124","B26","B34","B42","pr249a","pr248a","pr250a","pr242a","pr247a","183","218","181","198","239","220","PR160C","PR162C","PR164C","PR102C","101","PR109C","108","PR107C","106","PR143C","142","PR145C","PR151C","PR96C","95","PR147C","PR149C","PR105C","104","PR139C","PR141C"))
#View the tax_table
knitr::kable(head(tax_table(PR_Normal))) %>%
kableExtra::kable_styling("striped") %>%
kableExtra::scroll_box(width = "100%")
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| 1779c60c03d948a73687f440be973f68 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Flammeovirgaceae | Algivirga | Algivirga_pacifica |
| 25feebceed837c1dfa810cd8c0f26b41 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Cyclobacteriaceae | NA | NA |
| 0872b535e6e15b9cd031dae6a67d6257 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Amoebophilaceae | uncultured | uncultured_bacterium |
| 6590e42d5aeb18c338a836bc26028592 | d__Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | uncultured | uncultured | uncultured_Bacteroidetes |
| 0f22f07d74171ba4bbe6b3977217cbdc | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | uncultured_organism |
| d817bb1b166e13e13db6db8d56b43983 | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | NA |
#Now fix the labels
PR_Normal <- tax_fix(PR_Normal)
#View the relabeled tax_table
knitr::kable(head(tax_table(PR_Normal))) %>%
kableExtra::kable_styling("striped") %>%
kableExtra::scroll_box(width = "100%")
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| 1779c60c03d948a73687f440be973f68 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Flammeovirgaceae | Algivirga | Algivirga_pacifica |
| 25feebceed837c1dfa810cd8c0f26b41 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Cyclobacteriaceae | Cyclobacteriaceae Family | Cyclobacteriaceae Family |
| 0872b535e6e15b9cd031dae6a67d6257 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Amoebophilaceae | uncultured | uncultured_bacterium |
| 6590e42d5aeb18c338a836bc26028592 | d__Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | uncultured | uncultured | uncultured_Bacteroidetes |
| 0f22f07d74171ba4bbe6b3977217cbdc | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | uncultured_organism |
| d817bb1b166e13e13db6db8d56b43983 | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 Genus |
library(seecolor)
bcc_hex <- print_color(c("black", "#440154FF", "#450659FF","#460B5EFF","#472D7AFF","#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF","#2E6E8EFF", "#2D718EFF", "#2B748EFF","#2A778EFF", "#297B8EFF","#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF","#1FA287FF", "#21A585FF", "#23A983FF","#25AC82FF", "#29AF7FFF", "#2DB27DFF","#32B67AFF", "#37B878FF", "#3CBC74FF","#57C766FF", "#5EC962FF","#A2DA37FF","#DAE319FF", "#E4E419FF","#ECE51BFF", "#F5E61FFF","darkorchid1", "darkorchid2", "darkorchid3","#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF","#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#FF6699","#F68D45FF", "#FCA537FF","#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299","#300018"), type = "r")
##
## ------ c black #440154FF #450659FF #460B5EFF #472D7AFF #3B518BFF #3A548CFF #38598CFF #365C8DFF #34608DFF #33638DFF #31678EFF #306A8EFF #2E6E8EFF #2D718EFF #2B748EFF #2A778EFF #297B8EFF #1F958BFF #1F988BFF #1E9C89FF #1F9F88FF #1FA287FF #21A585FF #23A983FF #25AC82FF #29AF7FFF #2DB27DFF #32B67AFF #37B878FF #3CBC74FF #57C766FF #5EC962FF #A2DA37FF #DAE319FF #E4E419FF #ECE51BFF #F5E61FFF darkorchid1 darkorchid2 darkorchid3 #A21C9AFF #A62098FF #AB2394FF #AE2892FF #B22B8FFF #B6308BFF #BA3388FF #BE3885FF #C13B82FF #C53F7EFF #C8437BFF #CC4678FF #CE4B75FF #D14E72FF #D5536FFF #D7566CFF #DA5B69FF #DD5E66FF #E06363FF #FF6699 #F68D45FF #FCA537FF #F6E726FF #F4ED27FF #00489C #CCCCCC #999999 #A1C299 #300018 ------
## black
## #440154FF
## #450659FF
## #460B5EFF
## #472D7AFF
## #3B518BFF
## #3A548CFF
## #38598CFF
## #365C8DFF
## #34608DFF
## #33638DFF
## #31678EFF
## #306A8EFF
## #2E6E8EFF
## #2D718EFF
## #2B748EFF
## #2A778EFF
## #297B8EFF
## #1F958BFF
## #1F988BFF
## #1E9C89FF
## #1F9F88FF
## #1FA287FF
## #21A585FF
## #23A983FF
## #25AC82FF
## #29AF7FFF
## #2DB27DFF
## #32B67AFF
## #37B878FF
## #3CBC74FF
## #57C766FF
## #5EC962FF
## #A2DA37FF
## #DAE319FF
## #E4E419FF
## #ECE51BFF
## #F5E61FFF
## darkorchid1
## darkorchid2
## darkorchid3
## #A21C9AFF
## #A62098FF
## #AB2394FF
## #AE2892FF
## #B22B8FFF
## #B6308BFF
## #BA3388FF
## #BE3885FF
## #C13B82FF
## #C53F7EFF
## #C8437BFF
## #CC4678FF
## #CE4B75FF
## #D14E72FF
## #D5536FFF
## #D7566CFF
## #DA5B69FF
## #DD5E66FF
## #E06363FF
## #FF6699
## #F68D45FF
## #FCA537FF
## #F6E726FF
## #F4ED27FF
## #00489C
## #CCCCCC
## #999999
## #A1C299
## #300018
##Subsetting non-enriched samples to view composition in bar plot
PR_N <- PR_Normal %>% subset_samples(Sample_Type %in% c("N"))
PRSeqR_family <- PR_N %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Location,
Kingdom, Phylum, Class, Order, Family) %>%
filter(Abundance > 0.01) %>%
arrange(Class)
PRSeqR_family$Class_family <- paste(PRSeqR_family$Class, PRSeqR_family$Family, sep="_")
PRSeqR_family_bar <- ggplot(PRSeqR_family, aes(x = Location, y = Abundance, fill = Class_family)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +
guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") +
theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),
panel.grid.major = element_blank(), axis.text.x =element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),
axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),
axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),
legend.text = element_text(size = 7))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(PRSeqR_family_bar)
p <- ggplotly(PRSeqR_family_bar)
htmlwidgets::saveWidget(p, "PRSeqR_family_bar.html")
htmltools::tags$iframe(
src=file.path(getwd(), "PRSeqR_family_bar.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")
PRSeqR_Species <- PR_N %>%
tax_glom(taxrank = "Species") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Location,
Kingdom, Phylum, Class, Order, Family,Genus,Species) %>%
filter(Abundance > 0.05) %>%
arrange(Species)
PRSeqR_Species_bar <- ggplot(PRSeqR_Species, aes(x = Location, y = Abundance, fill = Species)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +
guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance (Species > 5%) \n") + xlab("Sample ID") +
theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),
panel.grid.major = element_blank(), axis.text.x =element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),
axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),
axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),
legend.text = element_text(size = 7))
p2 <- ggplotly(PRSeqR_Species_bar)
print(PRSeqR_Species_bar)
htmlwidgets::saveWidget(p, "PRSeqR_Species_bar.html")
htmltools::tags$iframe(
src=file.path(getwd(), "PRSeqR_Species_bar.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")
PRSeqR_family_Season <- PR_N %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Season,
Kingdom, Phylum, Class, Order, Family)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.01) %>%
arrange(Class)
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
## Creaitng the plot
PRSeqR_family_Season$Class_family <- paste(PRSeqR_family_Season$Class, PRSeqR_family_Season$Family, sep="_")
PRSeqR_Season_family_bar <- ggplot(PRSeqR_family_Season, aes(x = Season, y = Mean, fill = Class_family)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF","#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF","#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF","#2E6E8EFF", "#2D718EFF", "#2B748EFF","#2A778EFF", "#297B8EFF","#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF","#25AC82FF", "#29AF7FFF", "#2DB27DFF","#32B67AFF", "#37B878FF", "#3CBC74FF","#57C766FF", "#5EC962FF","#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF","darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF","#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#FF6699","#F68D45FF", "#FCA537FF","#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999","#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample ID") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_Season_family_bar)
p3 <- ggplotly(PRSeqR_Season_family_bar)
p
htmlwidgets::saveWidget(p, "Season_family.html")
htmltools::tags$iframe(
src=file.path(getwd(), "Season_family.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")
PRSeqR_family_Location <- PR_Normal %>%
tax_glom(taxrank = "Family") %>%
#Agglomerate at family level
transform_sample_counts(function(x) {x/sum(x)}) %>%
#Transform to rel.abundance
psmelt() %>%
# Melt to long format
group_by(Season, Location, Kingdom, Phylum, Class, Order, Family) %>%
dplyr::summarize(Mean = mean(Abundance, na.rm=TRUE)) %>% #Calculate average
filter(Mean > 0.01) %>%
#Filter arrange
arrange(Class)
## `summarise()` has grouped output by 'Season', 'Location', 'Kingdom', 'Phylum',
## 'Class', 'Order'. You can override using the `.groups` argument.
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_family_Location$Class_family <- paste(PRSeqR_family_Location$Class, PRSeqR_family_Location$Family, sep="_") #Creating plot
PRSeqR_Location_family_bar <- ggplot(PRSeqR_family_Location, aes(x = Location, y = Mean, fill = Class_family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7)) + facet_grid(. ~ Location, margins = TRUE, scale="free")
print(PRSeqR_Location_family_bar)
p4 <- ggplotly(PRSeqR_Location_family_bar)
htmlwidgets::saveWidget(p,"PRSeqR_Location_family_bar")
htmltools::tags$iframe(
src=file.path(getwd(),"PRSeqR_Location_family_bar"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")
PR_Normal <- PRphyseq %>% subset_samples(Sample_Type %in% c("N"))
PR_2022 <- PRphyseq %>% subset_samples(Year %in% c("2022"))
##Creating a subset of samples for only normal samples.This will exclude the enriched samples
PRSeqR_family_Year <- PRphyseq %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Year,Season, Kingdom, Phylum, Class, Order, Family)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.01) %>%
arrange(Class)
## `summarise()` has grouped output by 'Year', 'Season', 'Kingdom', 'Phylum',
## 'Class', 'Order'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_family_Year$Class_family <- paste(PRSeqR_family_Year$Class, PRSeqR_family_Year$Family, sep="_") #Creating plot
PRSeqR_Year_family_bar <- ggplot(PRSeqR_family_Year, aes(x = Year, y = Mean, fill = Class_family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Year") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7)) + facet_grid(. ~ Season, margins = TRUE, scale="free")
print(PRSeqR_Year_family_bar)
p5 <- ggplotly(PRSeqR_Year_family_bar)
Note that 2023 is only enriched samples so this year cannot be accurately depicted and 2022 is missing from the wet season.
*Now that we know there is an error in the 2022 data when attempting to generate the 2022 NMDS plot, we will be creating a family plot for the 2022 samples to determine what is too similar - this is what we think is preventing the NMDS plot from being generated.
PRSeqR_2022_Sample <- PR_2022 %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.00) %>%
arrange(Class)
## `summarise()` has grouped output by 'Location', 'Season', 'Kingdom', 'Phylum',
## 'Class', 'Order', 'Family', 'Genus'. You can override using the `.groups`
## argument.
##`summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_2022_Sample$Family <- paste(PRSeqR_2022_Sample$Class, PRSeqR_2022_Sample$Family, sep="_") #Creating plot
PRSeqR_2022_Sample_bar <- ggplot(PRSeqR_2022_Sample, aes(x = Location, y = Mean, fill = Family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_2022_Sample_bar)
p6 <- ggplotly(PRSeqR_2022_Sample_bar)
This is showing the samples averaged together, but there are multiple samples from each site. We are going to plot based on sample name to give us individual bars for each sample to better understand the similarities between each other.
PRSeqR_2022_Sample_ID <- PR_2022 %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Sample_ID,Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.00) %>%
arrange(Class)
## `summarise()` has grouped output by 'Sample_ID', 'Location', 'Season',
## 'Kingdom', 'Phylum', 'Class', 'Order', 'Family', 'Genus'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_2022_Sample_ID$Family <- paste(PRSeqR_2022_Sample_ID$Class, PRSeqR_2022_Sample_ID$Family, sep="_") #Creating plot
PRSeqR_2022_Sample_ID_bar <- ggplot(PRSeqR_2022_Sample_ID, aes(x = Sample_ID, y = Mean, fill = Family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample_ID") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_2022_Sample_ID_bar)
p7 <- ggplotly(PRSeqR_2022_Sample_ID_bar)
htmlwidgets::saveWidget(p,"PRSeqR_2022_Sample_ID_bar")
htmltools::tags$iframe(
src=file.path(getwd(),"PRSeqR_2022_Sample_ID_bar"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")
#PR_ENR <- PRphyseq %>% subset_samples(Sample_Type %in% c("E"))
##Creating a subset of samples for only enriched samples.This will exclude the normal samples
#PRSeqR_family_Location <- PR_ENR %>%
# tax_glom(taxrank = "Species") %>%
# transform_sample_counts(function(x) {x/sum(x)}) %>%
# psmelt() %>%
# group_by(Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%
# dplyr::summarize(Mean =
# mean(Abundance, na.rm=TRUE)) %>%
# filter(Mean > 0.00) %>%
# arrange(Class)
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_family_Location$Genus_species <- paste(PRSeqR_family_Location$Class, PRSeqR_family_Location$Species, sep="_") #Creating plot
## Warning: Unknown or uninitialised column: `Species`.
PRSeqR_Location_family_bar <- ggplot(PRSeqR_family_Location, aes(x = Location, y = Mean, fill = Genus_species)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Species > 0) \n") + xlab("Location") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_Location_family_bar)
p8 <- ggplotly(PRSeqR_Location_family_bar)
#(2/2) This step is needed to create an interactive plot on gtihub. After it has been pushed to github you can download the html doc and it will open as the interactive plot via your web browser.
htmlwidgets::saveWidget(p, "PRSeqR_Location_family_bar")
htmltools::tags$iframe(
src=file.path(getwd(), "PRSeqR_Location_family_bar"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0"
)
#Transform to relative abundance
PR_RA_Normal = transform_sample_counts(PR_Normal, function(x){x / sum(x)})
#View the otu_table to see the transformed relative abundances
View(otu_table(PR_RA_Normal))
#Choose colors
plot.colors <- c("steelblue2","purple4","darkorange","firebrick","springgreen4", "gold", "darkblue", "yellowgreen","turquoise4", "orange","indianred","darkslategrey", "lightblue","darkgreen","mediumaquamarine","gray48","mediumorchid1", "#5F7FC7","#DA5724", "#508578", "#CBD588","#CD9BCD","#AD6F3B", "#673770","#D14285", "#652926", "#C84248", "#8569D5", "#5E738F","#D1A33D","#8A7C64", "#599861","dodgerblue","darkmagenta", "forestgreen","steelblue1", "cyan","mediumorchid3", "cadetblue3", "yellow", "#00FF7F","#40E0D0","#E97451","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58", "#C4B454")
plot.shape <- c(21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21)
Subsetting from PR_physeq_RA_Normal, which includes only the non-enriched samples. This dataframe has been filtered so that ASVs that were 0 in all samples are excluded.
#PR_S003 <- PR_physeq_RA_Normal %>% subset_samples(Project %in% c("S003"))
#PR_S015 <- PR_physeq_RA_Normal %>% subset_samples(Project %in% c("S015"))
#PR_2021 <- PR_Normal %>% subset_samples(Year %in% c("2021"))
#PR_2022 <- PR_Normal %>% subset_samples(Year %in% c("2022"))
all.nmds.source.ord <- ordinate(
physeq = PR_RA_Normal,
method = "NMDS",
distance = "bray")
## Run 0 stress 0.1850681
## Run 1 stress 0.1850731
## ... Procrustes: rmse 0.0007928043 max resid 0.00578824
## ... Similar to previous best
## Run 2 stress 0.2139555
## Run 3 stress 0.2207744
## Run 4 stress 0.2022847
## Run 5 stress 0.2123585
## Run 6 stress 0.1850805
## ... Procrustes: rmse 0.001712935 max resid 0.01392611
## Run 7 stress 0.1871644
## Run 8 stress 0.1871644
## Run 9 stress 0.2024291
## Run 10 stress 0.2035301
## Run 11 stress 0.2080844
## Run 12 stress 0.1850833
## ... Procrustes: rmse 0.003884105 max resid 0.03295702
## Run 13 stress 0.1850062
## ... New best solution
## ... Procrustes: rmse 0.005042124 max resid 0.03300655
## Run 14 stress 0.2179661
## Run 15 stress 0.2297205
## Run 16 stress 0.1849911
## ... New best solution
## ... Procrustes: rmse 0.002372687 max resid 0.01513598
## Run 17 stress 0.1849877
## ... New best solution
## ... Procrustes: rmse 0.001707955 max resid 0.015439
## Run 18 stress 0.1850491
## ... Procrustes: rmse 0.005735748 max resid 0.05125967
## Run 19 stress 0.2237442
## Run 20 stress 0.1870003
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 18: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrmin
#Plot, color coding by particle type
all.nmds.Location <- plot_ordination(
physeq = PR_RA_Normal,
ordination = all.nmds.source.ord) +
scale_fill_manual(values = plot.colors, "Location")+
scale_shape_manual(values = plot.shape, name = "Location") +
geom_point(mapping = aes(fill = factor(Location), shape = Location, size = 5)) +
guides(size=FALSE) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
theme(plot.title = element_text(size = 18),
text = element_text(size = 18),
axis.title = element_text(size = 15),
panel.spacing = unit(1, "lines"),
panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_blank(),
legend.text = element_text(size = 15),
legend.title = element_text(size = 15),
legend.position = c(.99, .99),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
stat_ellipse(aes(group=Location))
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(all.nmds.Location)
#Plot, color coding by seaosn
all.nmds.Season <- plot_ordination(
physeq = PR_RA_Normal,
ordination = all.nmds.source.ord) +
scale_fill_manual(values = plot.colors, "Season") +
scale_shape_manual(values = c("D" = 21, "W" = 21), name = "Season") +
geom_point(mapping = aes(fill = factor(Season), shape = Season, size = 5)) +
guides(size=FALSE) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
theme(plot.title = element_text(size = 18),
text = element_text(size = 18),
axis.title = element_text(size = 15),
panel.spacing = unit(1, "lines"),
panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_blank(),
legend.text = element_text(size = 15),
legend.title = element_text(size = 15),
legend.position = c(.99, .99),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
stat_ellipse(aes(group=Season))
print(all.nmds.Season)
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ mm_Hg*DO_mg_L*DO_LDO*Specific_Conductance_uS_cm, data = filt_complete_bray_df)
#filt_complete
#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ Temperature_C*Location*pH*Salinity_ppt, data = filt_complete_bray_df)
#filt_complete
#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ Temperature_C*Location*DO_LDO*pH*Salinity_ppt, data = filt_complete_bray_df)
#filt_complete
#Convert output to dataframe
#filt.complete.pvalue <- as.data.frame(filt_complete)
#Create table of p-values
#knitr::kable(filt.complete.pvalue, row.names = NA) %>%
# kableExtra::kable_styling("striped",
# latex_options="scale_down") %>%
#kableExtra::scroll_box(width = "100%")
#Subset phyloseq object
#filt.mp.particle.RA <- subset_samples(MPfiltRA, sample_type == "Particle")
#sample_data(filt.mp.particle.RA)
#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_particle_bray <- phyloseq::distance(filt.mp.particle.RA, method = "bray")
#filt_particle_bray_df <- data.frame(sample_data(filt.mp.particle.RA))
#filt_all <- adonis2(filt_particle_bray ~ particle_type*effluent*week*polymer_type, data = filt_particle_bray_df)
#filt_all
#Convert output to dataframe
#filt.all.pvalue <- as.data.frame(filt_all)
#Create table of p-values
#knitr::kable(filt.all.pvalue, row.names = NA) %>%
#kableExtra::kable_styling("striped",
# latex_options="scale_down") %>%
#kableExtra::scroll_box(width = "100%")
#Run PERMANOVA with Bray-Curtis dissimilarity, only looking at particle type (glass vs. plastic)
#filt_particle_only_bray <- phyloseq::distance(filt.mp.particle.RA, method = "bray")
#filt_particle_only_bray_df <- data.frame(sample_data(filt.mp.particle.RA))
#filt_only_particle <- adonis2(filt_particle_only_bray ~ particle_type, data = filt_particle_only_bray_df)
#filt_only_particle
##Load libraries and read in your data
library(vegan)
library(ggplot2)
##Subset your data to only environmental data and abundance data.
##I already had a dataframe from the PR_RA_Normal for abundance data, so I just needed to do the same for the sample data (aka environmental data).
#df_env <- as.data.frame(sample_data(PR_RA_Normal))
#df_com <- as.data.frame(otu_table(PR_RA_Normal))
## Individual data frames created, now I am selecting the columns that will be used for this plot and renaming the data frame to com and env.
#com = df_com[,1:95]
#env = df_env[,15:17, 20, 21,24]
#convert com to a matrix and perform the NMDs ordination
#m_df_com = as.matrix(com)
#nmds code
#set.seed(123)
#nmds = metaMDS(m_df_com, distance = "bray")
#nmds
##Running envfit function with env dataframe
#en = envfit(nmds,env, permutations = 999, na.rm = TRUE)
library(phyloseq)
library(ggplot2)
library(microViz)
# install.packages("remotes")
#remotes::install_github("statdivlab/corncob")
#library(corncob)
#> microViz version 0.12.0 - Copyright (C) 2023 David Barnett
#> ! Website: https://david-barnett.github.io/microViz
#> ✔ Useful? For citation details, run: `citation("microViz")`
#> ✖ Silence? `suppressPackageStartupMessages(library(microViz))`
#knitr::opts_chunk$set(fig.width = 7, fig.height = 6)
#PR_RA_Normal <- corncob::PR_phylo
#PR_phylo
#> phyloseq-class experiment-level object
#> otu_table() OTU Table: [ 36349 taxa and 91 samples ]
#> sample_data() Sample Data: [ 91 samples by 15 sample variables ]
#> tax_table() Taxonomy Table: [ 36349 taxa by 7 taxonomic ranks ]
##Redundancy analysis is a constrained ordination method. It displays the microbial variation that can also be explained by selected constraint variables. Behind the scenes, a linear regression model is created for each microbial abundance variable (using the constraints as the explanatory variables) and a PCA is performed using the fitted values of the microbial abundances.
#PR_RA_Normal %>%
# ps_mutate(
# IBD = as.numeric(ibd == "ibd"),
# Female = as.numeric(gender == "female"),
# ) %>%
# tax_transform("clr", rank = "Genus") %>%
# ord_calc(
#constraints = c("IBD", "Female", "Abx."),
# method = "RDA", # Note: you can specify RDA explicitly, and it is good practice to do so, but microViz can guess automatically that you want an RDA here (helpful if you don't remember the name?)
# scale_cc = FALSE # doesn't make a difference
# ) %>%
# ord_plot(
# colour = "DiseaseState", size = 2, alpha = 0.5, shape = "active",
# plot_taxa = 1:8
# )